Tarea 1: Foundations

CC6104: Statistical Thinking

Integrantes :

Cuerpo Docente:

Índice:

  1. Objetivo
  2. Instrucciones
  3. Referencias
  4. Segunda Parte: Elaboración de Código

Objetivo

a la primera tarea del curso Statistical Thinking. Esta tarea tiene como objetivo evaluar los contenidos teóricos de la primera parte del curso, los cuales se enfocan principalmente en análisis exploratorio de datos y conceptos introductorios de probabilidades. Si aún no han visto las clases, se recomienda visitar los enlaces de las referencias.

La tarea consta de ejercicios prácticos con el fín de introducirlos a la programación en R enfocada en el análisis estadístico de datos.

Instrucciones:

Referencias:

Slides de las clases:

Videos de las clases:

En la siguiente sección deberá resolver cada uno de los experimentos computacionales a través de la programación en R. Para esto se le aconseja que cree funciones en R, ya que le facilitará la ejecución de gran parte de lo solicitado.

Pregunta 1: Visualización de Datos

Para esta pregunta usted deberá trabajar en base al conjunto de datos hearth_database.csv, el cual esta compuesto por las siguientes variables:

En base al dataset propuesto realice un análisis exploratorio de los datos (EDA). Para su análisis enfoquen el desarrollo en las siguientes tareas:

Respuesta

Pregunta 2: Teorema Central del Limite

Pruebe el teorema central del limite aplicando un muestreo de la media en las distribuciones Gamma, Normal y una a su elección. Grafique los resultados obtenidos y señale aproximadamente el numero de muestreos necesarios para obtener el resultado esperado, pruebe esto con las siguientes cantidades de muestreo \(\{10,100,1000,5000\}\). ¿El efecto ocurre con el mismo número de muestreo para todas las distribuciones?.

Respuesta

# Definición de variables o estructuras necesarias para el muestreo.
# Primero definimos la función plot_gauss, que nos permitirá graficar rápidamente.
plot_gauss <- function(x, name) {
  avg <- mean(x)
  std <- sd(x)
  x_axis <- seq(min(x), max(x), length = length(x))
  y_axis <- dnorm(x_axis, mean = avg, sd = std)
  hist(x, main = sprintf("Gráfico de la distribución %s", name))
  lines(x_axis, y_axis, type = "l", col = "red")
}

# Representamos las distribuciones Normal, Gamma y de Poisson, de parámetros N(1.5, 0.75), G(2, 1.5) y P(5).
avg <- 1.5
std <- 0.75
shape <- 2
scale <- 1.5
lambda <- 5

# Inicializamos los vectores que acumularán las medias para cada muestreo.
iter <- 1000
normal_vector <- vector(length = iter)
gamma_vector <- vector(length = iter)
poisson_vector <- vector(length = iter)
# Graficamos las 3 distribuciones para distintos tamaños de la muestra
sample_size <- 10

# Realizar el muestreo de las distribuciones.
for(i in 1:iter) {
  new_vector <- rnorm(n = sample_size, mean = avg, sd = std)
  normal_vector[i] <- mean(new_vector)
  
  new_vector <- rgamma(n = sample_size, shape = shape, scale = scale)
  gamma_vector[i] <- mean(new_vector)

  new_vector <- rpois(n = sample_size, lambda = lambda)
  poisson_vector[i] <- mean(new_vector)
}

# Código para gráficos
par(mfrow = c(1, 3))
plot_gauss(normal_vector, sprintf("Normal con tamaño %s", sample_size))
plot_gauss(gamma_vector, sprintf("Gamma con tamaño %s", sample_size))
plot_gauss(poisson_vector, sprintf("de Poisson con tamaño %s", sample_size))

# Graficamos las 3 distribuciones para distintos tamaños de la muestra
sample_size <- 100

# Realizar el muestreo de las distribuciones.
for(i in 1:iter) {
  new_vector <- rnorm(n = sample_size, mean = avg, sd = std)
  normal_vector[i] <- mean(new_vector)
  
  new_vector <- rgamma(n = sample_size, shape = shape, scale = scale)
  gamma_vector[i] <- mean(new_vector)

  new_vector <- rpois(n = sample_size, lambda = lambda)
  poisson_vector[i] <- mean(new_vector)
}

# Código para gráficos
par(mfrow = c(1, 3))
plot_gauss(normal_vector, sprintf("Normal con tamaño %s", sample_size))
plot_gauss(gamma_vector, sprintf("Gamma con tamaño %s", sample_size))
plot_gauss(poisson_vector, sprintf("de Poisson con tamaño %s", sample_size))

# Graficamos las 3 distribuciones para distintos tamaños de la muestra
sample_size <- 1000

# Realizar el muestreo de las distribuciones.
for(i in 1:iter) {
  new_vector <- rnorm(n = sample_size, mean = avg, sd = std)
  normal_vector[i] <- mean(new_vector)
  
  new_vector <- rgamma(n = sample_size, shape = shape, scale = scale)
  gamma_vector[i] <- mean(new_vector)

  new_vector <- rpois(n = sample_size, lambda = lambda)
  poisson_vector[i] <- mean(new_vector)
}

# Código para gráficos
par(mfrow = c(1, 3))
plot_gauss(normal_vector, sprintf("Normal con tamaño %s", sample_size))
plot_gauss(gamma_vector, sprintf("Gamma con tamaño %s", sample_size))
plot_gauss(poisson_vector, sprintf("de Poisson con tamaño %s", sample_size))

# Graficamos las 3 distribuciones para distintos tamaños de la muestra
sample_size <- 5000

# Realizar el muestreo de las distribuciones.
for(i in 1:iter) {
  new_vector <- rnorm(n = sample_size, mean = avg, sd = std)
  normal_vector[i] <- mean(new_vector)
  
  new_vector <- rgamma(n = sample_size, shape = shape, scale = scale)
  gamma_vector[i] <- mean(new_vector)

  new_vector <- rpois(n = sample_size, lambda = lambda)
  poisson_vector[i] <- mean(new_vector)
}

# Código para gráficos
par(mfrow = c(1, 3))
plot_gauss(normal_vector, sprintf("Normal con tamaño %s", sample_size))
plot_gauss(gamma_vector, sprintf("Gamma con tamaño %s", sample_size))
plot_gauss(poisson_vector, sprintf("de Poisson con tamaño %s", sample_size))

Notemos que para la distribución de Poisson, los datos convergen a una distribución Normal al realizar muestreos de tamaño 100, mientras que las otras dos aún no muestran ese comportamiento. Es al tamaño de muestra de 1000 que las 3 distribuciones se comportan como uno esperaría teóricamente.

Pregunta 3: Ley de los Grandes Numeros.

Lanzamiento de monedas

Realice el experimento de lanzar una moneda cargada 1000 veces y observe el comportamiento que tiene la probabilidad de salir cara. Para realizar el experimento considere que la moneda tiene una probabilidad de \(5/8\) de salir cara. Grafique el experimento para las secuencias de intentos que van desde 1 a 1000, señalando el valor en que converge la probabilidad de salir cara.

Respuesta

# Simular lanzamientos
for (lanzamientos in 1:1000) {
  
}

# Gráfico de la convergencia

El problema de Monty Hall

Remontándonos en la televisión del año 1963, en USA existía un programa de concursos donde los participantes debían escoger entre 3 puertas para ganar un premio soñado. El problema del concurso era que solo detrás de 1 puerta estaba el premio mayor, mientras que detrás de las otras dos habían cabras como “premio”.

Una de las particularidades de este concurso, es que cuando el participante escogía una puerta, el animador del show abría una de las puertas que no fue escogida por el participante (Obviamente la puerta abierta por el animador no contenía el premio). Tras abrir la puerta, el animador consultaba al participante si su elección era definitiva, o si deseaba cambiar la puerta escogida por la otra puerta cerrada. Un vídeo que puede ayudar a comprender mejor este problema en el siguiente link.

Imagine que usted es participante del concurso y desea calcular la probabilidad de ganar el gran premio si cambia de puerta en el momento que el animador se lo ofrece. Utilizando listas/arrays/vectores simule las puertas del concurso, dejando aleatoriamente el premio en alguna posición del array. Hecho esto, genere un numero de forma aleatoria para escoger una de las puerta (posiciones de la estructura), para luego ver si cambiando de posición tendrá mayores posibilidades de ganar el premio. Genere N veces el experimento y grafique cada una de las iteraciones, tal como se hizo en el ejercicio de las monedas.

Respuesta:

# Creamos una función que simule el juego
montyhall <- function(cambiar = TRUE){
  Puertas <- sample(1:3,3)             #Puertas donde la posición que tiene el 3 es el premio
  posicion <- sample(1:3,1)            #Elección del participante.
  
  return(Eleccion) # Retornamos la elección, esta puede que tenga el premio o no
}

# Función que simula N juegos
n_juegos <- function(n = 10 ,cambiar_puerta = TRUE){
  
}

Pregunta 4: ¿Independencia?

Ustedes disponen de los dados D1 y D2, los cuales no están cargados y son utilizados para comprobar que \(\mathbb{P}(AB)=\mathbb{P}(A)\mathbb{P}(B)\) cuando el evento A es independiente del B. Para estudiar la independencia considere que los eventos A y B se definen de la siguiente manera; sea A el evento dado por los valores obtenidos en el lanzamiento del dado D1, este está compuesto por \(A=\{D1=1,D1=2,D1=6\}\). Por otro lado, el evento B viene dado por los valores obtenidos con el dado D2, el que está conformado por \(B=\{D2=1,D2=2,D2=3,D2=4\}\). Con esto, tendremos un \(\mathbb{P}(A)=1/2\), \(\mathbb{P}(𝐵)=2/3\) y \(\mathbb{P}(AB)=1/3\). Compruebe de forma gráfica que al realizar 1000 lanzamientos (u otro valor grande que usted desea probar) se visualiza que \(\mathbb{P}(AB)=\mathbb{P}(A)\mathbb{P}(B)\).

Hecho lo anterior, compruebe el comportamiento de un segundo grupo de eventos, dados por el lanzamiento de solo el dado D1. Donde, los eventos para D1 quedan definidos como: \(A =\{D1=1,D1=2,D1=6\}\) y \(B=\{D1=1,D1=2,D1=3\}\). ¿Se observa independencia en este experimento?.

Se le aconseja que para simular los lanzamientos de dados utilice la función sample() para generar valores aleatorios entre 1 y 6. Compruebe los números generados por la función con los casos favorables de cada uno de los eventos a ser estudiados.

Respuesta:

n_lan <- 6000000

# Primer grupo de eventos
dice1 <- sample(1:6, n_lan, replace = TRUE)  # Dado 1
dice2 <- sample(1:6, n_lan, replace = TRUE)  # Dado 2
dice_count1 <- table(dice1)
dice_count2 <- table(dice2)

favorable_a <- c(1, 2, 6)  # {1, 2, 6}
favorable_b <- c(1, 2, 3, 4)  # {1, 2, 3, 4}
favorable_ayb <- c(1, 2)  # {1, 2}

count_a <- sum(dice_count1[favorable_a])  # |A|
count_b <- sum(dice_count2[favorable_b])  # |B|

# Ahora contamos los casos donde A y B se cumplen al mismo tiempo
count_ayb <- 0  # |A∩B|
for (i in 1:n_lan) {
  if (dice1[i] %in% favorable_a && dice2[i] %in% favorable_b) {
    count_ayb <- count_ayb + 1
  }
}

# Calculamos las probabilidades
prob_a <- count_a / n_lan  # P(A)
prob_b <- count_b / n_lan  # P(B)
prob_ayb <- count_ayb / n_lan  # P(A∩B)

print(sprintf("P(A) = %s", prob_a))
## [1] "P(A) = 0.4998125"
print(sprintf("P(B) = %s", prob_b))
## [1] "P(B) = 0.666653833333333"
print(sprintf("P(A)*P(B) = %s", prob_a*prob_b))
## [1] "P(A)*P(B) = 0.333201919072917"
print(sprintf("P(A∩B) = %s", prob_ayb))
## [1] "P(A∩B) = 0.333163833333333"

Ahora visualizamos las probabilidades

par(mfrow = c(1, 2))
y_ax <- c(0, 0, 1, 1)
x_ax <- c(0, 1, 1, 0)

# Graficamos P(A)
plot(1, ann = FALSE, type = "n", axes = FALSE, xlim = c(0, 1), ylim = c(0, 1), main = "P(A)")
x_p_a <- c(0, prob_a, prob_a, 0)
x_p_a_c <- c(prob_a, 1, 1, prob_a)

polygon(x_p_a, y_ax, col = "blue")
text(mean(x_p_a), mean(y_ax), expression(P(A)), col = "white", cex = 1.5)

polygon(x_p_a_c, y_ax, col = "red")
text(mean(x_p_a_c), mean(y_ax), expression(P(A^c)), col = "white", cex = 1.5)


# Graficamos P(B)
plot(2, ann = FALSE, type = "n", axes = FALSE, xlim = c(0, 1), ylim = c(0, 1), main = "P(B)")
y_p_b <- c(prob_b, prob_b, 1, 1)
y_p_b_c <- c(0, 0, prob_b, prob_b)

polygon(x_ax, y_p_b, col = "blue")
text(mean(x_ax), mean(y_p_b), expression(P(B)), col = "white", cex = 1.5)

polygon(x_ax, y_p_b_c, col = "red")
text(mean(x_ax), mean(y_p_b_c), expression(P(B^c)), col = "white", cex = 1.5)

Ahora grafiquemos \(\mathbb{P}(AB)\)

plot(1, ann = FALSE, type = "n", axes = FALSE, xlim = c(0, 1), ylim = c(0, 1), main = "P(A∩B)")

polygon(c(0, prob_a, prob_a, 0), c(prob_b, prob_b, 1, 1), col = "blue")
text(prob_a / 2, 1 - prob_b / 4, expression(P(A * "∩" * B)), col = "white", cex = 1.5)

polygon(c(0, prob_a, prob_a, 0), c(0, 0, prob_b, prob_b), col = "purple")
text(prob_a / 2, prob_b / 2, expression(P(A * "∩" * B^c)), col = "white", cex = 1.5)

polygon(c(prob_a, 1, 1, prob_a), c(prob_b, prob_b, 1, 1), col = "purple")
text(1 - prob_a / 2, 1 - prob_b / 4, expression(P(A^c * "∩" * B)), col = "white", cex = 1.5)

polygon(c(prob_a, 1, 1, prob_a), c(0, 0, prob_b, prob_b), col = "red")
text(1 - prob_a / 2, prob_b / 2, expression(P(A^c * "∩" * B^c)), col = "white", cex = 1.5)

Notemos que el área del cuadrado de la esquina superior izquierda, que representa el área de \(\mathbb{P}(A \cap B)\)

Pregunta 5: La Ruina del Jugador

Un amigo ludópata suyo le comenta que el truco de jugar en el casino esta en no parar de apostar y apostando lo mínimo posible. Ya que así, tienes mas probabilidades de ganar el gran pozo que acumula el juego. Usted sabiendo la condición de su amigo, decide no creer en su conjetura y decide probar esto a través de un experimento.

Para realizar el experimento usted decide asumir las siguientes declaraciones, bajo sus observaciones:

En el primer experimento deberá obtener la evolución de los fondos hasta que el jugador se queda sin fondos para jugar. Puede ser útil seguir la lógica de una moneda cargada para realizar esto. Pruebe esto con una apuesta igual a 5, 25 y 50 graficando los resultados. Comente los resultados obtenidos.

Para la segunda parte de este experimento, con las funciones ya creadas, proyecte 5000 apuestas y obtenga la probabilidad a la que converge el experimento empezando con una apuesta de: 5, 25 y 50. Para este experimento considerar como éxito (1) cuando su función ruina supera los 200 y considere lo contrario cuando su función perdida es menor o igual a 0.

Respuesta

# Función para obtener el desarrollo de las apuestas
ruina <- function(money = 100, bet = 5) {
  vec_money <- c(money)
  while (0 < money && money < 200) {
    game <- as.integer(runif(1, 1, 19))
    sign <- if (game <= 8) 1 else -1
    money <- money + sign * bet
    vec_money <- append(vec_money, money)
  }
  return(vec_money)
}

plot(ruina(), type="l", col="blue", xlab="N° de juegos", ylab="Fondos", main="Evolución de los fondos (apuesta = 5)")

plot(ruina(bet = 25), type="l", col="blue", xlab="N° de juegos", ylab="Fondos", main="Evolución de los fondos (apuesta = 25)")

plot(ruina(bet = 50), type="l", col="blue", xlab="N° de juegos", ylab="Fondos", main="Evolución de los fondos (apuesta = 50)")

En general se ve que la tendencia es perder, solo que al apostar el mínimo se pierde más lento ya que al ser una variación más pequeña, se desestabiliza (es decir, se va a los extremos) de a poco, no así al apostar 50, que es posible perder en dos juegos.

n_bet <- 5000
bet5_vector <- vector(length = n_bet)
bet20_vector <- vector(length = n_bet)
bet50_vector <- vector(length = n_bet)
for (i in 1:n_bet) {
  final5_money <- tail(ruina(bet = 5), 1)
  final20_money <- tail(ruina(bet = 20), 1)
  final50_money <- tail(ruina(bet = 50), 1)
  win5 <- if (final5_money <= 0) 0 else 1
  win20 <- if (final20_money <= 0) 0 else 1
  win50 <- if (final50_money <= 0) 0 else 1
  bet5_vector[i] <- win5
  bet20_vector[i] <- win20
  bet50_vector[i] <- win50
}
t5 <- table(bet5_vector)
t20 <- table(bet20_vector)
t50 <- table(bet50_vector)
print(sprintf("Apuesta = 5 -> E = %d, F = %d, Prob = %f", t5["1"], t5["0"], t5["1"]/(t5["1"] + t5["0"])))
## [1] "Apuesta = 5 -> E = 62, F = 4938, Prob = 0.012400"
hist(bet5_vector, main = sprintf("Éxitos y fracasos para apuesta = 5"))

print(sprintf("Apuesta = 20 -> E = %d, F = %d, Prob = %f", t20["1"], t20["0"], t20["1"]/(t20["1"] + t20["0"])))
## [1] "Apuesta = 20 -> E = 1220, F = 3780, Prob = 0.244000"
hist(bet20_vector, main = sprintf("Éxitos y fracasos para apuesta = 20"))

print(sprintf("Apuesta = 50 -> E = %d, F = %d, Prob = %f", t50["1"], t50["0"], t50["1"]/(t50["1"] + t50["0"])))
## [1] "Apuesta = 50 -> E = 1938, F = 3062, Prob = 0.387600"
hist(bet50_vector, main = sprintf("Éxitos y fracasos para apuesta = 50"))

 

A work by CC6104